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Abstract: We propose a new method for the numerical solution of back- 

^.^ ward stochastic differential equations (BSDEs) which finds its roots in 

«^ Fourier analysis. The method consists of an Euler time discretization of the 

JT^ BSDE with certain conditional expectations expressed in terms of Fourier 

r^ transforms and computed using the fast Fourier transform (FFT). The 

(-H problem of error control is addressed, we consider the extension of the 

.4.^ method to reflected BSDEs, and some numerical examples are considered 

C^ from finance demonstrating the performance of the method. 
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1. Introduction 



—il Backward stochastic differential equations (BSDEs) have been a topic of interest 

f^ since the early work of Bismut [4] and the results of Pardoux and Peng [25] on 

C^ their well-posedness. A BSDE is an equation of the form 

^ — I 

.> Yt^C+ f{s,Y,,Zs)ds- ZldW, (1.1) 

k> Jt Jt 

j^ defined on a complete filtered probability space (fi, V,F, {^t}t^[o.T]) where W 

is a standard n-dimensional Brownian motion, the terminal condition ^ £ K*^ is 
a square integrable J-r-measurable random variable and the driver / : [0, T] x 
MJ^ is a functional. It is known from [25] that there exists a 



X 



pfcxr 



unique adapted square integrable backward process Y taking values in M*^ and 
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a unique predictable process Z with values in M"^*^ satisfying equation (1.1) 
under Lipschitz and integrability conditions on the driver /. 

Many works have extended this existence and uniqueness result. Antonelli [1] 
introduced forward-backward stochastic differential equations (FBSDEs). Mao 
[23], Lepeltier and San Martin [20] and Kobylanski [19] among others treat 
non-Lipschitz cases. Also, the theory of BSDEs has found various applications 
particularly in finance and in the study of partial differential equations (PDEs). 
From Pardoux and Peng [26] (see also El Karoui, Peng and Quenez [14, Section 
4.1]) we have that if the Cauchy problem to the one-dimensional diffusion PDE 

iw + k&+ /(^'"' if) = 0' (i'^) G [0'^) X ^ (1.2) 

[u(T,x) = g{x) 

has a unique solution u e C^'^ then the solution (Y, Z) for the BSDE (1.1) with 
terminal condition ^ ~ 5(M^t) admits the representation 

Yt = u{t,Wt) (1.3) 

Zt = £it^Wt). (1.4) 

Conversely, the solution of the PDE (1.2) can be interpreted in terms of the 
solution of the BSDE (1.1). General formulations of the nonlinear Feynman- 
Kac formula for FBSDEs, quasilinear parabolic PDEs, and viscosity solutions 
have been studied extensively. 

Deriving an explicit solution to a nontrivial (F)BSDE is possible only in very 
few situations, such as Yong [31], Hyndman [18] and Richter [29]. Thus, numer- 
ical methods for BSDEs have been studied extensively. Numerical methods for 
(F)BSDEs can be classified into three main groups: PDE based methods, spa- 
tial discretization based methods, and Monte-Carlo based methods. PDE based 
methods, which started with the finite difference approach of Douglas, Ma and 
Protter [12], consider a numerical resolution of the nonlinear parabolic PDE 
related to the (F)BSDE. The two other methods rely on a time discretization 
of the (F)BSDE. Spatial discretization based methods (see Chevance [10], Bally 
and Pages [2], Delarue and Menozzi [11] or Peng and Xu [27] among others) 
use a deterministic space grid. On the other hand, the space discretization is 
random in Monte-Carlo based methods (for instance, Bouchard and Touzi [7], 
Gobet, Lemor and Warin [17], and Bender and Denk [3]). 

In this paper, we propose an alternative spatial discretization method for 
BSDEs and illustrate its implementation in the one-dimensional case. To the 
best of our knowledge, the most efhcient approach in this simple case is the 
binomial method of Peng and Xu [27] which has connections with the theoretical 
work of Briand, Delyon and Memin [8] and Ma et al. [22]. However, our method 
avoids a notable drawback of the binomial method: the contraction of the space 
grid leading to the approximation of the Wiener process by means of scaled 
random walks. Indeed, we use a fixed equidistant space grid, thus allowing an 
exact simulation of the Wiener process at time nodes. The EFT algorithm. 
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which plays a key role in our method, helps in producing an efficient algorithm. 
As in Carr and Madan [9] and Lord et al. [21] in the context of option pricing 
under Levy processes, we employ the FFT algorithm to compute quadratures. 
The presence of dynamic programming through the Euler scheme is a major 
similarity between our method and Lord et al. [21]. 

This paper is structured as follows. Section 2 reviews Euler time discretization 
schemes for BSDEs which are used in Section 3 to develop the convolution 
method. Section 4 presents a detailed error analysis of the convolution method. 
Some extensions of the method are presented in Section 5, numerical results for 
examples from finance are included in Section 6, and Section 7 concludes. 

2. Time discretization of BSDEs 

The convolution method developed in this paper, as any spatial discretization 
based method, requires the availability of a time discretization scheme for BS- 
DEs. In this section, we present the Euler time discretization schemes that are 
widely used in numerical methods for BSDEs. Alternatives to the Euler schemes 
can be found in the 0— schemes of Zhao, Chen and Peng [34] or the penalization 
scheme proposed by Peng and Xu [27] inspired by a method used by El Karoui 
et al. [16] to prove well-posedness of reflected BSDEs. Convergence of the Euler 
schemes are considered by Zhang [32, 33] and Bouchard and Touzi [7]. 

For simplicity of notation we shall suppose all processes are one-dimensional 
{k = n = 1). Further, we make the following assumption to ensure existence 
and uniqueness of the BSDE (1.1). 

Assumption 2.1. We suppose a Markovian terminal condition with 

C = g{WT) (2.1) 

where g : R ^ R is real function satisfying the square integrability condition 

E [e] = E [giWrf] < ^. (2.2) 

In addition, both the terminal condition g and the driver f verify the Lipschitz 
condition 

\g{x) - g{x)\ + |/(i, y, z) - f{i, y, z)\ < C {\x - x\ + \y - y\ + \z - z\) (2.3) 

for some constant C > 0, Vx, x, y, y, z,z ^ M, and Vt, t > 0. 

Consider the time mesh tt = {0 = tg < ii < ... < in = T} on the time 
interval [0, T] with a number n e N of time steps. The Euler scheme leads to 
the following time discretization 

Yt: = Y,^^ + f{t,,Y,-^, Zl)A, ~ ZlAW,. (2.4) 

for the BSDE of equation (1.1) where A^ — i^+i — ti is the time step and 
AWi — Wti^i — Wti is the Brownian increment. Moreover, Y^ and Z^. stand for 
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the approximate solution at time node ti. Taking the conditional expectation of 
both sides of equation (2.4) yields 

Y,: = E [y,:^^ I J-,,] + fiU, Yl, Zl)A,. 

If we first multiply both sides of equation (2.4) by the Brownian increment AWi 
and then take the conditional expectation we find 



zi 



1 

a; 



E 



Y^,AWAJ^t 



ti 



i\-rt. 



The backward algorithm for numerical solution of BSDEs is then defined by 



ZJ" = 0, Fj'^ = C 



Zl 



ilE 



Y^^AWATti 



< i <n 



(2.5) 



^*: = E 



Yt:j:Fu +/(i„i;:,z-)A„o<z< 



and is known as the implicit Euler scheme since the value of the approximate for- 
ward process Y^ appears on both sides of the last equation. Zhang [33] proposes 
choosing ^^ so that the quadratic approximation error is of first order 



Air II 2 



where |7r| is the maximal time step 



le-n' =o{\n\) 



(2.6) 



\7t\ := max \ti+i - U 

0<i<n 



(2.7) 



However, in the context of BSDEs with a Markovian terminal condition, we set 

C = ^ = giWr) . 

In order to avoid solving a non-linear system of equations to recover the 
approximate forward process values, one may consider an alternative scheme 
that is explicit in those values 



zi = 0, Y,i - r 

y/:aw.\^u 



K = s-E 



Yu=^ 



, <i <n 



(2.8) 



r,-^^ + /(i„y,-^^,z-)A,|j-,, 



, < i < n 



which we call the explicit Euler scheme I. Another explicit scheme consists of 
replacing the conditional expectation of the driver in the explicit Euler scheme 
I by the driver evaluated at the conditional expectations of the arguments. This 
procedure leads to 





= 0, 


Yt: 


. = ^" 








zi 


^ilE 


Y/:^^AW.\J^u 


, <i < n 


Yu 


= E 


YtU^u 


+ fiU,E 


Yt:j^u_ 



(2.9) 



,Zl)A,,0<i<n 
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which we call the explicit Euler scheme II. The approximate {Y, Z) processes 
then takes the form 

Yt^ = r,:, Z[ = zi for t e [UM+i)- (2.10) 

on the entire time interval. 

The global discretization error E.^ is defined as 



Ei := max E 

0<i<ri-l 



sup \Yt-Yl^ 
«e[ti,t.+i] 



4=0 



+ 1 



Z,-Zl? ds 



(2.11) 



for any version of the Euler scheme. The following theorem, from Zhang [33] or 
Bouchard and Touzi [7], gives an error order for the implicit Euler scheme. 

Theorem 2.1. Under the setting of Assumption 2.1, the implicit Euler scheme 
yields a first order quadratic error 

K=0{\7r\). (2.12) 

Due to the Lipschitz nature of the driver /, it can be proved by induction 
that the quadratic error between the implicit scheme and the explicit schemes 
is of first order. As a consequence, the explicit schemes also have a first order 
quadratic error as noted by Bouchard, Elie and Touzi [6, Remark 2.1.1]. 

3. Convolution method 

In this section, we detail the convolution method for the conditional expectations 
involved in the Euler scheme. We also present the discretization of these quadra- 
tures and their relationship to the discrete Fourier transform (DFT) which can 
be efficiently computed using the fast Fourier transform (FFT). 

3.1. Convolution on the explicit Euler scheme II 

The starting point of the convolution method for BSDEs is an explicit Euler 
scheme. If we consider the explicit Euler scheme II of equation (2.9) an approxi- 
mate solution of the BSDE (1.1) at mesh time ti consists of real- valued functions 
Mi, iii, and iii defined by the backward recursions 

Mj(x) = Ui{x) + Aj/(ii, Ui{x),Ut{x)) (3.1) 

Mx) = T^E [u,+iiWt,^JAW,\Wt, = x] 

"a"/ (y - ^)'^^+iiy)Hy - x)dy (3.2) 

u,[x)^^[u,+^{Wu^,)\Wu =x] 

Ut+i{y)h{y - x)dy (3.3) 
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for i = 0, 1, ..., n— 1 and m„(x) — g{x). Note that Ui represents the approximate 
Y process and iii stands for the approximate Z process at mesh time ti while 
Ui is an intermediate quantity. The notation u used in equations (3.1)-(3.3) is 
not to be confused with the solution of the PDE (1.2) as we do not employ the 
representation (1.3)-(1.4) in this paper. The function h is the density function 
of Wf.^, conditional on the value of Wt. 

h{x)^{27TA,)--^exp(^-^y (3.4) 

If a method for calculating the integrals of equations (3.2) and (3.3) is avail- 
able, then the sequence {ui{Wti) , UiiWtJ) for i = 0,1, 2, ...,n — 1 is an approxi- 
mation to the BSDE solution of equations (1.3) and (1.4) on the interval [0,T]. 
The stationarity and independence of Brownian increments allow us, as in Lord 
et al. [21], to express the functions Ui and iii in equations (3.2) and (3.3) as con- 
volutions. These convolutions suggest the using Fourier transforms and hence 
the computation of the integrals via discrete Fourier transforms. 

Recall the Fourier transform of an integrable real function rj is the function 
?7 : M — > C is defined as 

/oo 
e-'''^f]{x)dx (3.5) 

-oo 

where i = \/— T is the imaginary unit. The inverse Fourier transform recovers 
the function r/ from its Fourier transform fj through the relation 

l{x) := ^-\fi]{x) - — / e'--f^{v)dv. (3.6) 

In many cases of interest the function 77 is not integrable. For any real function 
?7 : M — > M define the dampened function 77" as 

ri°^{x) = e-"^77(x). (3.7) 

where a € M is a dampening parameter chosen so that 77" is integrable. 
Taking the Fourier transform of itf in equation (3.3) gives 



:SK]{v)^ I e-'-'^e-^M u,+i{y)h{y - x)dydx 

y J —00 

e""'" / u'^,{y)e"^y--'>h{y-x)dydx 



— CX3 



= 5K+i](^)5[e-"^/i(-^)](^) (3.8) 

using the convolution theorem. Moreover, making the change of variable x = —2;, 

/oo pOO 

e-'"''e-°"'h{-z)dz = / e'(''-'")^/7(a;)dx 
00 J —00 

= </>(t/-ia) (3.9) 
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where 0(i^) — exp (— ^AjZ/^) is the characteristic function of the density h. 

The equahty of equation (3.9) is well-defined since \(j){i^ ~ '^ct)\ < oo for any 
a G M. Nonetheless, the structure of the terminal condition g (and more gener- 
ally, the preceding approximation ut^i) will have a major impact in the choice of 
the dampening parameter a. Indeed, combining equations (3.8) and (3.9) gives 

i?[<](i^)=5[<+i](^.)0(^.-i«) (3.10) 

and hence the parameter a must be chosen so that the dampened functions uf , 
i = 0,l,...,n, are integrable and admit Fourier transforms. 

Similarly, the Fourier transform of iif in equation (3.2) is given by 

1 

A" 



^[<](^) = -— ;?K+i](^)J[ze'"^/i(-z)](^) 



-^5[<+i](^)|;5[e-"^M-^)](^) 



-— 5K+i]HT^'^(^-i«) 



= {a + iiymu2+i]{i^)cj){j^ - ia) (3.11) 

using the differentiation properties of the Fourier transform. 

From equations (3.10) and (3.11), we recover the functions Ui and iii by 
taking the inverse Fourier transform and adjusting for the dampening factor 

u^{x) = e"^;?-i[;?K+i](z.)0(z.-ia)](a;) (3.12) 

Mx) = e""5-i[(a + M5K+i](z/)0(z/-ia)](a;). (3.13) 

Equations (3.1), (3.12), and (3.13), evaluated at a; = Wt^, define a convolution 
method for the approximate solution of the BSDE (1.1) based on the explicit 
Euler scheme II. 

3.2. Convolution on the explicit Euler scheme I 

An alternative characterization of the approximate solution of the BSDE (1.1) 
is obtained if one considers the explicit Euler scheme I of equation (2.8). In this 
case, the approximate solution {Y, Z) consists of functions Vi and Vi at mesh 
time ti which take the form 

v,{x) ^'E[v,+i{Wt^^,)\Wt^ ^ x] 

Vi+i{y)h{y - x)dy (3.14) 

J — oo 

where 

Ui+i{x) = Vi+i{x) + Aif{ti,Vi+i{x),Vi{x)), (3.15) 

v,{x) = ^E [v,+i{Wu^,)AW,\Wu = x] 
1 [°° 

"" TT / iy- ^)^^+iiy)Hy-^)dy, (3.i6) 



C.B. Hyndman and P. Oyono Ngou/Convolution method for BSDEs 8 

for i — 0, 1,...,T7- — 1 and Vn{x) — g(x). 

Following the steps of the previous characterization equations (3.14) and 
(3.16) lead naturally to 

v,{x) = e"^5-i [S[vf+^]{v)^{v - la)] [x) (3.17) 

v,{x) = e"^;?-i [{a + iv):S[vf+^]{i^)<P{v ~ la)] {x). (3.18) 

where both vf and vf for i = 0, l,...,n— 1 along with the dampened terminal 
condition are assumed to be integrable so that they admit Fourier transforms. 
Equations (3.15), (3.17), and (3.18) define a convolution method for the approx- 
imate solution of the BSDE (1.1) based on the explicit Euler scheme I. 

3.3. Numerical implementation 

From equations (3.12) and (3.13) or equations (3.17) and (3.18) one notices that 
computing the approximate solutions (u^, Ui) and (w^, iii) at mesh time t,i reduces 
to computing a function : M ^ M depending on two functions f/' : C — > C and 
77 : M — ^ M in the following manner 

e{x) = — I e"'''^{v)xl;{v)dy (3.19) 

if we drop the dampening factor e"^ . 

This integral is numerically computed by discretizing the Fourier space with 
a uniform grid of A^ + l points {I'ilflo ^^ ^^'^ interval [— f , f ] of length L, where 
N is even, such that 

v., = fo + iAi/ (3.20) 

with vq — —^ and /S.v — ^. Hence, for any a; G M 

L 

e{x) « — /' e"'''^{v)tlj{v)dv 
27tJ_l 

« ^ E ^''''''?"(^^)'^(^^) (3-21) 

where the integral is approximated using lower Riemann sums and 

/•OO /"OO 

^{ly,) = / e-''"'^'n°'{x)dx = / e-'^'''e-°^r?(a;)da;. (3.22) 



00 



This last integral is also computed using an uniform grid of A'^+l points {xj}^q 
on the restricted interval [a;o,a;Ar] centered at W^ = such that 

Xj = --+jAx (3.23) 

where Ax = ^^ is chosen so that the Nyquist relation LI = 2'kN is satisfied. 
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The discretization of the integral in equation (3.22) leads to an expression 
involving the discrete Fourier transform (DFT). The DFT is a numerical pro- 
cedure that transforms a set of real or complex numbers {xj},Sq^ into another 
set {xj}fSQ^ through the relation 



i.:=2)[a;], = ^^e-^'=^x, (3.24) 

for k = 0,1,. ..,N — 1. The inverse DFT performs a reciprocal operation by 
computing the set of numbers {xj}jSq using the numbers {xj}jSq as 

w-i 
Xk := T)-'[x]k = J2 e'^'^'^xj (3.25) 

j=o 

for fc = 0, l,...,iV- 1. 

As we intend to use the DFT to compute (3.21) and (3.22) we assume that 
the following conditions are satisfied. 

Assumption 3.1. The generic dampened function rj" and its derivative -4j— 
have the same values at the boundaries of the restricted domain [xq,xm] 

77"(xo) = r{xN), (3.26) 

%ix,) . %ix.). (3.27) 

We approximate the integral of equation (3.22) by first restricting the inte- 
gration interval to [xq, xn] and then applying a composite quadrature rule with 
weights {ifilflo- Consequently, 

^{vi) « f " e~''"'^'n°'{x)dx (3.28) 

N 

« AxJ2w,e-'^^''^r,^{x,) 

N 

= Ax ■ 6-'^°'^' Y^ Wje-'^'^e-'^''"^''jf{xj) 

3=0 

2n 



e--o..j, [{{~iyw,f^-ix,)}fSo']^ (3.29) 



for i = 0, 1, ...,7V — 1 using Assumption 3.1 since N is even and e^"'"'^^ = —1 
with 

Wj — Wj + Sn-j,nwn (3.30) 

where Sij stands for the Kronecker's delta. At this point, many quadrature 
rules are available. For example, one may use the composite trapezoidal rule 
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with weights of the form 

w,^l-]^{5o,. + 5N.i) ,i = 0,l,...,iV (3.31) 

leading to Wj = 1. Higher order composite quadrature rules will improve accu- 
racy in presence of a smooth driver /. 

A similar approach can be found in Lord et al. [21] who enhance the discrete 
Fourier transform with a composite trapezoidal quadrature rule to compute 
this last integral. However, [21] omits Assumption 3.1 making the error analysis 
quite tedious and leading to considerable numerical errors, especially around 
the boundaries of the restricted domain. In addition, we use a fixed space grid 
whereas [21] shift the space grid through time steps. The major difference, how- 
ever, is that the scheme in [21] solves the Snell envelope and does not seek a 
numerical solution for BSDEs. 

Remark 3.1. From equation (3.28), the restriction of the real line to the in- 
terval [xq,xn] also solves the problem of integrability of the generic dampened 
function rj" that we pointed out on equation (3.10). Indeed, the restriction is 
essentially equivalent to considering the truncated function rj^l^^^.^^-^. Hence, 
Tj" needs to be integrable only on the restricted domain. 

The values of the functions 9 are computed at the grid points {xkJi^^Q by 
combining equations (3.21) and (3.29) 

N-l 

e{xu) « Y. e''^^^''V(^,)e-'^°''^D [{(-l)^^.ry«(x,)},^o'], 

N-l 



j=o 



{-\f^- 






(3.32) 



Since we use the DFT, the underlying trigonometric (and hence periodic) inter- 
polation allows us to set 

e{xN)^e{xo). (3.33) 

In applications we shall consider functions 77" that do not satisfy Assump- 
tion 3.1. To address this problem we slightly modify the function rj by adding a 
linear function to obtain a modified dampened function 77^ ^ defined as 

r^%,{x)^e-'^-{r^{x)+px + K). (3.34) 

The following lemma gives the appropriate choice for the dampening parameter 
a € M, and the coefficients /3 S M and k £ M. 

Lemma 3.1. Suppose the real function rj G C^[a,b] is differentiable with 
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and let ri^ ^ be the modified, dampened function defined in equation (3.34)- Then 
a ^ ^^oJp^ (3.35) 



e'"''iv{b) + I3b) - e-°"'{T]{a) + /3a) 



e~' 
solve the system of nonlinear equations 



(3.36) 



V^pja) = rj'^pjb) (3.37) 



for anyf3i{-^{a),-l^{b)}.In addition, if 



>max(|g(6)|,|g(a)|) (3.39) 



then a G M and k G 



Proof. Equation (3.37) gives (3.36) in a straightforward manner using basic 
algebra. Since rj is difFerentiable rj^ ^ is also differentiable and 



dx 



{x) = -ar^%,{x)+e-'^^l^^{x)+l3 



and combining both equations (3.37) and (3.38) leads to (3.35). Clearly, if the 
inequality (3.39) holds then both (g^(&) + /3) and (g^(a) + /3) are strictly positive 
and a € M. D 

Remark 3.2. When implementing the method, the values of the derivative ^^ 
at xq and x^ can be approximated by forward or backward finite differences. 

Remark 3.3. A positive constant, e > 0, which represents the minimal slope 
allowed in the linear transform j3x + k is used to ensure equation (3.39) is 
enforced. Set 

/? = e + max('||^(x^)|,||^(xo)|V (3.40) 

Whenever -(^{xn) = 7^{^o)> one can set a = k = Q and 

f]{xN)~1l{xo) 



xn - xo 



(3.41) 



Under the transformation of equation (3.34), the computation of our ap- 
proximate solution is not significantly more complex. Since the function is a 
(dampened) conditional expectation, properties of the conditional expectation 
allows the adjustments as shown by the following Theorem. 
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Theorem 3.1. Let 77 : [a, 6] — > M be an integrable Junction and define rj^a '■ 
[a, b]^^ as 

such thatrj'^ ^ is the modified, dampened function ofrj defined by equation (3.34)- 
Then the function : [a, 6] — ?> M of equation (3.19) admits the alternative repre- 
sentation 

^i^)^^J e"'^^Jiy)ij;{iy)du^H{x,a,P,K) (3.42) 

where 

Hix, a, f3, n) = {^""^^ '^ "^^"^ = ^^ + ''''^^^'' " ^"^ (3.43) 

^ ^ \e-""(/3a; + K) if ^{ly) ^ (j){i^ - ia). ^ ' 

Proof. First, let ip{i') = {a + iiy)(/){iy — ia). By definition, we know that 
e{x) ^ ^ I e'''''^{u)^{u)dv 

ZTT 



e 



-00 

'Q.X 



A,, 



A,, 



- (E [{v{Wu+,) + PWu+, + At) ^W,\Wu =x\- /3A, 



-E [r]p^,{Wu^,)AW,\Wu =x]- e-""/3 



A, 
^ ^/^e'^-7^,(z.)^(z.)di.-e-"-/3. 

Similarly, if 'ip{i') = (piv — ia), we have 

e{x) ^ ^ f e""'^{v)^{v)dv 

271' J ^00 



e-^''^[7^{Wu^,)\Wu=x\ 

e-"^E [vp.,^{Wu^,)\Wu =x]- e-"^(/3x + n) 

^ r e'^'^^^(v)i,{v)dv-e-^'^{Px + n). 



n 



Theorem 3.1 shows that the computational formula of equation (3.32) can 
be applied to the dampened transform function ryS ^ which satisfies Assumption 
3.1 after an appropriate choice of the coefficients a, (3 and n using Lemma 3.1. 
One recovers the function 9 by subtracting the function H . 

The implementation of the convolution method gives the approximation val- 
ues {uik}k=o ^'^'i {'^ifclfeLo to the approximate solutions Ui and iii for i = 
0, 1, 2, ..., n — 1. An approximate solution to the BSDE consists of a (linear) in- 
terpolation of a Brownian path through the node values {uik}^^^ and {uifej^o 
for J = 0, 1, 2, ..., rt — 1. We next consider the problem of computational error. 
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4. Error Analysis 

The convolution method induces two main types of error. In addition to the time 
discretization error E^^ discussed in Section 2, there is a space discretization error 
and we focus on this last error term. We limit our analysis to the explicit Euler 
scheme II since equivalent results are obtained for the explicit Euler scheme I 
using the same techniques. 

Throughout the section, {uifejj^o' {"«fc}fcLo ^'^'^ {'"ifel^o denote the numer- 
ical solution of equation (3.32) obtained from the convolution method at time 
mesh ii, i = 0, 1, ..., n — 1. The convolution method induces a space discretiza- 
tion error when approximating the values of Ui{xk) and Ui{xk) by un- and Uik 
respectively. We will particularly describe the local behavior of this error term. 
We define it as 

Eik := \ui{xk) - Uik\ + \uiixk) ~Uik\- (4.1) 

The following lemma describes the DFT accuracy in approximating the Fourier 
coefficients and proves useful in the derivation of a space discretization error 
bound. We skip the proof since the results are well known (see Plato [28, The- 
orem 3.4, p. 140] and Vretblad [30, Theorem 4.4, p. 85]). 

Lemma 4.1. Suppose the integrable function u : [—212] ^ ^ with 
u(— 2) = "(9) o,dm,its the Fourier series expansion 



u{x) = 



CO 

E 

A:— — 00 



Cke ' 



^k^ 



X e 



I I 

2' 2 



(4.2) 



and {xk}k=o '^'"'^ *^^ nodes of the equidistant grid of [— |, |] such that 
Xk = —^ + fcA where TV G N is even and A = N~^l. 
IfueC^ then 

Ck.^^{-l)''~^V[{i-iyuix,)}fs,']^ + OiA^) 

for fc = 0, 1, ..., A'^ — 1 and 

\ck\<Ck-^ , 

for k e Z\{0} and some constant C > depending on j^ . Consequently. 

^-1 



uix) 



E 



CkC 



ik^x 



-©(A), Vxe 



k=- 



l I 

2' 2 



(4.3) 
(4.4) 

(4.5) 



The next theorem gives an error bound for the space discretization error 
under smoothness conditions on the BSDE coefficients / and g. 

Theorem 4.1. Suppose f E C^'^'^ and g E C^ . Then for any i = 0, 1, ...,n — 1 
and k = 0,1, ..., N , the convolution method applied on the truncated interval 
[—i, I] yields a (local) discretization error of the form 



E^k = x{xk) + O (Ax) + O (e-^^"'') 



(4.6) 
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where the extrapolation error x satisfies 



\x{xk)\<C\ / h{y)dv\ (4.7) 

for some positive constants C,K > Q depending on the driver f , the function g, 
and the terminal time T when using the trapezoidal quadrature rule. 

Proof. Suppose the solution u^+i at time t^+i is known. Since / e C^'^'^ and 
g G C^, it is easily shown that Ui+i e C^. Also, we know from Zhang [33] and 
Bouchard and Touzi [7] that Y^ ^ = Ui+i{Wt-^^) is square integrable so that 
Ui+i is square integrable (with respect to the Gaussian density). 

In the light of Theorem 3.1, we can limit ourselves to the case where 

f l\ fl\ , 9-Uj+i / l\ du,+i fl 

so that a ~ P = n ^ {). Let {cfc}^_^ be the Fourier coefficients of Ui+i on 



Ui{xk) = / Ui+i{y)h{y - Xk)dy + i Ui+i{y)h{y - Xk)dy 



where 



Ui+i{xk + y)h{y)dy + / Ui+i{xk + y)h{y)dy 

\y\<h •'\v\>i 



/ Ui+i{xk + y)h{y)dy = E Ui+i{xk + AWn-i)lM.\l-^,^]i^^n-i] 



for some constant i^ > by successively applying Cauchy-Schwartz and Cher- 
noff inequalities since the solution u^+i is square integrable. Hence 

u,{xk) = / Too{xk + yMy)dy + o(e-''^^'"'") 
J\y\<h ^ ' 

+ {ui+i{xk + y) -Too{xk + y))h{y)dy (4.8) 

Ay\<^ 

where T^oix) — J2'kL-oo cte^^^^ for cc G M. So that, on one hand, we have 
Tcaixk + y)h{y)dy 

= Y. c^e'^'^^-cj^fj^)^ f T^{xk + y)h{y)dy 



\y\<i 
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= E cje'^''^^'=<^ (j^) - f T^ixk + y)h{y)dy + 0{Ax) 
,_« V ' / J\y\>i 

J~ 2 

(by Lemma 4.1), 
„■__« \ / 

J~ 2 

(by boundedness of Too and ChernofF inequality), 

JV-l 

= (-1)'= ^ 0(z/j)(-l)^-^c^_«e'T^^''= + 0{Ax) + O (e-^^'"''') 

JV-l 

+ ©(Aa;) + O {e-^^^''^^\ (by Lemma 4.1), 
= u,fc+C'(Aa;)+o(e-^'^' "'''). (4.9) 

Assuming x^ > 0, without loss of generality, define xo as 



Xo(a;fe) = / {ui+i{xk + y)~Too{xk+y))h{y)dy 

J\y\<i 

{ui+i{xk + y)~ u.i+i{xk +y -I)) h{y)dy 

since Too is periodic and Too(x) — Ui+x{x) on the interval [^5 '^2]- Equation 
(4.8) becomes 

u,[xk) = u,k + Xo{xk) + 0{Ax) + O (e-^^'"'') (4.10) 

and we note, by the continuity of Wi+i, that 



\Xo{xk)\ < Co / hiy)dy (4.11) 



for some positive constant Co > 0. 
Similar techniques show that 



u^ixk) = u,k + Xi{xk) + O (Ax) + O (e-^^'"''') (4.12) 



C.B. Hyndman and P. Oyono Ngou/Convolution method for BSDEs 



16 



where K >Q and, letting v{y) = ^,+1(2;^ + y) - T^{xk + y), 

Xi{xk) ^ A~^ / yv{y)h{y)dy 
J\y\<i 

J\y\<i y 



= a; 



y 



\y\<i 



dv 
dx 



(y) 



dx^ 



iOy h{y)dy 



\y\<h 



y'^ -Q^{y)Hy)dy 



= A: 



Since ^ 

OX 



\y\< 






(for some ^ G 
(by symmetry) 
y) ] h{y)dy. 



l_ l_ 
2' 2 



is the Fourier expansion of ^^^^ , we have 




y^hiy)dy 



(4.13) 



by the continuity (and boundedness) of "'1'^'" and Cauchy-Schwartz inequahty, 
for some constant Ci > 0. 

The Lipschitz property of the driver / completes the proof from the relations 
in equations (4.9), (4.11), (4.12) and (4.13). D 

Theorem 4.1 decomposes the spatial discretization error in three parts: the 
truncation error, the discretization error and the extrapolation error. Most PDE 
based and spatial discretization based methods for BSDEs fail in giving a bound 
for the error due to truncation. The error analysis shows that the truncation 
error 0(e^^' ) has a spectral convergence of index 2 when applying the con- 
volution method. Also, the discretization error O (Ax), of first order, is similar 
to PDE based methods such as Douglas, Ma and Protter [12] or Milstein and 
Tretyakov [24]. 

The extrapolation error x is specific to the convolution method implemented 
using the DFT. Equation (4.7) shows that errors appear and may accumulate 
around the boundaries of the truncated domain. Nonetheless, the extrapolation 
error is mainly time related through the density h and can be confined at the 
boundaries for fine time discretizations as shown in the following corollary. 



Corrolary 4.1. Under the conditions of Theorem 4-1, 

lim x{xk) = 



(4.14) 



for any Xk G (-|, |). 
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Proof, li Xk = then equation (4.7) gives |x(0)| < and the result holds. If 
Xk ^ and x^ e (^1:1): then 



lim / h{y)dy 



lim / h{y)dy 



lim / h{y)dy 



-\xk\ 



S{y)dy 



(where 6 is the Dirac delta function), 
= 

since ^ [^ — \xk\ , |] . Equation (4.7) then leads to the result. 



D 



5. Extensions 

Various simple extensions can be made of the convolution method. One of the 
most important is reflected BSDEs. We also consider the convolution method 
under arithmetic Brownian motion. These cases have interesting applications in 
mathematical finance, especially for option pricing. 

5.1. Reflected BSDEs 

Explicit Euler schemes have been constructed for reflected BSDEs with contin- 
uous barrier which make it possible to apply the convolution method to such 
BSDEs. Consider the solution (Y, Z, A) of the system 



-dYt = fit, Yt, Zt,)dt + dAt - ZtdWt 
Yt>Bt ,dAt>0,yte[0,T] 
[J^{Yt - Bt)dAt = , Ft = 9{Xt) 



where the lower barrier is a deterministic function B : [0, T] x 
and the Brownian motion and 

Bt = B{t,Wt). 

This RBSDE is associated to the following obstacle problem 

u{t,x) > B{t,x), {t,x) e [o,r] X M 



(5.1) 



of time 



(5.2) 



(5.3) 



u{T,x) =g(a;), x £ 
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as established by El Karoui et al. [16] . An adaptation of the explicit Euler scheme 
I provides the numerical solution to the reflected BSDE through the equations 



(Zl^O.Ytl^C 



Zl = StE 



y^^Aw-.IJ-,^ 



AAl = E y^^ +/(f„y,7^^,Z-)A,|J-,, -BiU,Wu) 



(5.4) 



'-ti 



E 



r- +/(i.,r,-^^,z-)A,|j-i, 



^K 



where for any number a; € M, a;^ = max(0, —a;). 

The problems of time discretization of RBSDEs and their convergence were 
treated in Bouchard and Chassagneux [5] for the implicit Euler scheme. Peng 
and Xu [27] proposed an equivalent scheme with a discrete filtration and proved 
its convergence under a binomial method. The scheme is easily solved with a 
convolution method by first noticing that the approximate solution (u^, u^, Ai^j) 
at mesh time ii, where , Vi is the approximate reflection process, can be written 
as 



where 



Vi{x) = 
Vi+l{x) = 

Avt{x) = 



Vi+i{y)h{y - x)dy + Avi{x) 



Vi+i{x) + Aif{U,Vi+i{x),v^{x)) and 
{y - x)vi+i{y)h{y - x)dy 

Vt+i{y)Hy - x)dy - B{U,x) 



(5.5) 

(5.6) 

(5.7) 

(5.8) 



for i = 0,l,...,n— 1 and Vn{x) — g{x). The computation of Vi and the integral 
part of the approximate solution Vi is identical to the non-reflected case. 

One can also naturally build an alternative scheme from the explicit Euler 
scheme II. The approximate solution u^, the approximate gradient Ui and the 
approximate reflection Aui at mesh time ti then take the form 



Ui{x) 



where 



Ui{x) + Aif{ti,Ui{x),Ui{x)) + Au^{x) 



Ui{x) = ^ (y - x)ui+i{y)h{y - x)dy 

/oo 
Ui+iiy)h{y - x)dy 
CO 

Aui{x) = {ui{x)+ Aif[ti,Ui[x),Ui{x))- B{U,x)Y 



(5.9) 

(5.10) 

(5.11) 
(5.12) 



for i = 0, 1, ..., n — 1 and Un{x) = g{x). 
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5.2. Arithmetic Brownian motion and Levy processes 

We can extend the convolution method to consider an arithmetic Brownian 
motion 

Xt = xo + iit + aWt (5.13) 

as the forward process and solve the associated Cauchy problem for the advection- 
diffusion equation 

[iE + i^Tr. + 5^'& + m^:",'^^) = , (i,x) e [0,T) X M 
\u{T,x) =g{x), xe M 

to which an obstacle can be added when in presence of a reflected BSDE. The for- 
ward process increments are stationary, independent, and normally distributed 
with density 

h{x) = ^ exp (J^-J^^\ (5.15) 

and characteristic function 

<^(y)^e^^^'^'''-^y''^). (5.16) 

The development of the convolution method in this case also leads to transforms 
identical to equation (3.19) with 

i}{v)^(i)(y-\a) (5.17) 

when computing the approximate solutions Wj or the intermediate solutions Ui 
and 

^{v)=a{a^\v)(i)(y-\a) (5.18) 

when computing iii and Vi. In our implementation, the ii and v are actually 
estimates for (t|j but the scheme can easily be modified so as to estimate the 
derivative |^ directly. 
The equivalence 

1 Z""" . . 

271' 7-00 

= ^/ e'''''^Mi^i^)dy-H[x,a,P,K) (5.19) 

of Theorem 3.1 still holds with 

H(x,a, p, Kj = < _ 

le °''-^l3a ii xp[i/) — a [a + ih')4){h' — ia). 

(5.20) 
Also, the convolution method can be used to compute conditional expecta- 
tions under general Levy processes as in Lord et al. [21]. Indeed, the indepen- 
dence of increments and the availability of the characteristic function are the 
only requirements to apply the method. The convolution method may also serve 
as a numerical method for partial differential integral equations (PIDE) under 
a Levy process. 
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6. Application to option pricing 

Through option pricing problems, we will particularly treat the case of BSDEs 
with non-linear and non-smooth drivers that present the lowest rate of con- 
vergence. An introduction to financial applications of BSDEs, particularly to 
imperfect markets and American option problems, can be found in El Karoui 
and Quenez [15], El Karoui, Pardoux and Quenez [13] or El Karoui, Peng and 
Quenez [14]. 

For the market model consisting of a single risky asset (or stock) S with the 
dynamics 

St = e^* (6.1) 

where the process X represents the stock return. We first consider an European 
call option with maturity T = 1 and strike price K under a lending rate of 
r = 0.01 and a borrowing rate R. The return process is an arithmetic Brownian 
motion 

Xt=Xo+L-6-^aAt + <jWt (6.2) 

such that the stock has an initial value of So = e^° = 100, an expected return 
rate of // = 0.05, a dividend rate of 5 and a volatility of cr = 0.2. 

The European call option price then follows a BSDE with the return process 
X as the forward process, the driver 

/(t, y, z) = -ry - (^) z + (i? - r) (y - ^) " (6.3) 

and the terminal function 

g{x) = (e- - Kf (6.4) 

under the given imperfect market conditions. The American call option solves 
a reflected BSDE with the barrier function 

B(t, x) = g{x) = (e' - K)+ , {t, x) e [0, T] x R. (6.5) 

When the borrowing rate equals the lending rate R — r — 0.01 and S — 0, the 
European and the American call options have the same price. Figure 1 shows the 
structure of the absolute log error on European option prices and deltas where 
the true values are computed using the Black-Scholes formula. As expected the 
errors are amplified at the boundaries of the truncated domain, but also for 
around-the-money options, to a lesser extent, due to the non-smoothness of the 
terminal function g. In addition, out-of-the-money options have smaller absolute 
errors compared to in-the-money options. 

The Black-Scholes formula gives call option prices of 4.6101, 8.4333 and 
14.1929 at strike prices K = 110, 100 and 90 respectively when the other vari- 
ables are kept unchanged. The true values for the option deltas are 0.3720, 
0.5596 and 0.7507 when the strike price is K = 110, 100 and 90 respectively. 
Table 1 gives the relative error in percentage on price estimates for the Amer- 
ican call options with both convolution schemes using different time steps and 
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Fig 1. Absolute error of the convolution method (scheme II) for American call option prices 
and deltas (non-dividend paying stock and no market frictions) 




[xo,xm] = Xo + [-5,5], Af = 2i2, n = 1000, e = 5, K = So = 100, R = r = 0.01, <5 = 0. 



the indicated strike prices. Table 2 gives the estimates for the American option 
deltas obtained from the approximate gradient by 



Delta 






(6.6) 



when using the explicit Euler scheme II. The option prices and option deltas 
are valued on the restricted domain [xq, xm] = Xq + [—5, 5] with N = 2^^ grid 
points and a minimal slope of e = 5 (see Remark 3.3). 

The results of Table 1 and 2 show the accuracy of the convolution method 
on a RBSDE with a smooth linear driver. Indeed, the relative error percentages 
remain low (less than 0.3%) for the estimated option prices and deltas. However, 
out-the-money option estimates seem to display the largest relative errors. 

For a borrowing rate of _R = 0.03 and a dividend rate of (5 = 0, Table 3 shows 
the estimates for the American call prices when the option is at the money 
Sq = K — 100 and r = 0.01. Moreover, the convolution methods return an 
option delta of 0.5987 when applied with n ~ 2000 time steps. In this case, the 
Black-Scholes formula does not apply but since the American and the European 
call options have the same price, one can notice that the simulated paths never 
hit the barrier in Figure 2. Paths are simulated using the solution from the 
convolution method applied on scheme II on the restricted domain [xq , xn] = 
Xq + [—5, 5] with N = 2^^ grid points, n = 1000 time steps and e = 5. We used 
n = 1000 time steps to simulate the stock price [St). 
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Table 1 
Relative errors (%) for American call option prices (non-dividend paying stock and no 

market frictions) 





K (Strike) 


n=500 


n=1000 


n=2000 


n=5000 




110 


0.0456 


0.0217 


0.0108 


0.0043 


Convolution (Scheme I) 


100 


0.0178 


0.0095 


0.0047 


0.0024 




90 


0.0049 


0.0028 


0.0014 


0.0007 




110 


0.0087 


0.0239 


0.0022 


0.0001 


Convolution (Scheme II) 


100 


0.0059 


0.0024 


0.0012 


0.0007 




90 


0.0028 


0.0014 


0.0007 


0.0004 



[xo.xjv] =X„ + [-5,5], N ■- 



e = 5, R = r = 0.01, <S = 0. 



Table 2 
Relative errors (%) for American call option deltas (non-dividend paying stock and no 

market frictions) 



K (Strike) 


90 


100 


110 


Convolution (Scheme I) 


0.0133 


0.0010 


0.2414 


Convolution (Scheme II) 


0.0133 


0.0010 


0.2414 



[xo,xn] =Xo + [-5,5], N ■- 



5, R- 



:0.01, (5 = 0. 



Table 3 
ATM American call option prices (non-dividend paying stock under imperfect market 

conditions) 

n (number of time steps) 500 1000 2000 5000 

Convolution (Scheme I) 9.4127 9.4131 9.4132 9.4133 
Convolution (Scheme II) 9.4132 9.4133 9.4133 9.4134 



[xo.xjv] =Xo + [-5,5], N: 



5, R = 0.03, r = 0.01, K = So = 100, <5 = 0. 
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Fig 2. (ATM) American call option sample paths (non-dividend paying stock under imperfect 
market conditions) 




S 40 



a 20 




g 0.5 

■i, 



t (time) 




Pathi 
Path 2 









Paths 















0.5 1 

t (time) 

[xo.xjv] =X„ + [-5,5], N = 212, „^ ;^oOO, e = b, K = So = 100, iJ = 0.03, r = 0.01, <5 = 0. 



The option price can be calculated using a Monte-Carlo method such as 
the forward scheme of Bender and Denk [3]. However, in the context of uni- 
dimensional BSDEs, Monte-Carlo methods will generally be less efficient than 
space discretization methods. As an illustration, the convolution method on 
both explicit Euler schemes runs in approximately 4 seconds when pricing the 
options of Table 3 with n — 1000 time steps. On the other hand, the forward 
scheme runs in 18 seconds with only n = 20 time steps. We used the 7 first 
power functions as basis functions and 100, 000 paths to generate the Monte- 
Carlo estimates. The Picard iterations are stopped whenever the difference in 
two consecutive prices is less than 10~^ for a maximum number of 10 iterations. 
Fifty (50) independent valuations with the Monte-Carlo method give a 95% 
confidence interval of [9.3972, 9.4222] which includes all estimates of Table 3. 
Hence, the convolution method gives satisfactory results even for coarse time 
discretization. 

Table 4 provides price estimates for out-of-the-money and in-the-money op- 
tions and Table 5 gives estimates for option deltas on a non-dividend-paying 
stock under imperfect market conditions. Both tables compare the estimates 
obtained with the convolution method and those obtained with the binomial 
method of Peng and Xu [27]. The convolution method and the binomial method 
give similar prices and delta values for all options which confirms the convolu- 
tion method accuracy even for non-smooth drivers. Nonetheless, the binomial 
method is faster (less that tenth of a second for 1000 time steps) than the con- 
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Table 4 
American call option prices (non-dividend paying stock under imperfect market conditions) 



n (number of time steps) 


K (Strike) 


500 


1000 


2000 


5000 


Convolution (Scheme I) 


110 
90 


5.2924 
15.4289 


5.2929 
15.4291 


5.2931 
15.4292 


5.2933 
15.4292 


Convolution (Scheme II) 


110 
90 


5.2932 
15.4290 


5.2933 
15.4291 


5.2933 
15.4291 


5.2934 
15.4292 


Binomial Method 


110 
90 


5.2918 
15.4313 


5.2945 
15.4263 


5.2936 
15.4298 


5.2937 
15.4292 



[xo.xjv] = Xq + [-5,5], N = 2^2, e = 5, R = 0.03, r = 0.01, 5 = 0. 

Table 5 
American call option deltas (non-dividend paying stock under imperfect market conditions) 



K (Strike) 


90 


100 


110 


Convolution (Scheme I) 


0.7814 


0.5987 


0.4104 


Convolution (Scheme II) 


0.7814 


0.5987 


0.4104 


Binomial Method 


0.7813 


0.5987 


0.4104 



[xo, xjv] = Xo + [-5, 5], AT = 212, e = 5, n = 200O, R = 0.03, r = 0.01, 5 = 0. 



volution method for the same num.ber of time step when computing the BSDE 
initial values. 

If we introduce a dividend rate oi S = 0.035 under imperfect market condi- 
tions {R = 0.03 and r = 0.01), then the American and the European call option 
prices differ and the Black-Scholes formula does not apply. The convolution 
method estimates the (at-the-money) American call option price at 7.5610 and 
the European call option price at 7.4712. We use scheme II with the restricted 
domain [xo,xn] = ATq + [—5,5], N ~ 2^^ grid points, n = 2000 time steps 
and a minimal slope of e = 5. Figure 3 shows the typical sample paths for the 
American option where the reflecting process At (hedging cost) is now non-zero 
for in-the-money path indicating a difference in price with the European call 
option. 

Figure 4 displays the option price and delta surfaces. The regularity of these 
surfaces indicates that the convolution method is efficient in handling non- 
smoothness in the terminal condition g but also in the driver /. 

7. Conclusion 

We considered explicit Euler time discretizations to develop a new spatial dis- 
cretization method for backward stochastic differential equations. This new 
method is based on the expression of conditional expectations as convolutions 
and the FFT algorithm is used to compute those integrals and as such we refer 
to it as the convolution method. Since the FFT algorithm is more suitable for 
periodic functions we introduced a transform in order to treat BSDEs with non- 
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Fig 3. (ATM) American call option sample paths (dividend-paying stock under imperfect 
market conditions) 




0.5 
t (time) 

[xo,X]v] = Xo + [-5,5], N = 2^2, „, = 1000, e = b, K = So = 100, R = 0.03, r = 0.01, 

5 = 0.035. 



Fig 4. (ATM) American call option price and delta surfaces 
(dividend paying stock with market frictions) 



a) Price surface 




S (st()cl< price) 







b) Delta surface 



t (time) 








S (stocl< price) 



t (time) 



[xo, Xm] = ^0 + [-5, 5], Af = 2^2, n = 1000, e = b, K = So = 100, R = 0.03, r = 0.01, 

5 = 0.035. 
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periodic terminal conditions. We then provided a (local) error analysis which 
indicates that the use of the FFT performs an extrapolation that induces a 
non-negligible error term meaning that transform does not completely solve 
the problem of non-periodicity. We extended the method to consider reflected 
backward stochastic differential equations and equations driven by arithmetic 
Brownian motion which are important classes of BSDEs in finance. 

Numerical experiments show that the convolution method is accurate and 
handles non-linearity and non-smoothness in the BSDE coefficients. The ad- 
dition of a technique to suppress the extrapolation error is an interesting im- 
provement to the method and shall be presented in a future paper. Finally, 
further generalizations of the method to forward-backward stochastic differen- 
tial equations and efficient implementations for multidimensional problems are 
important areas of future research. 
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